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Supercritical Flow Past a Symmetrical Bicircular Arc Airfoil 


Abstract 

Li and Holt (1981) developed a numerical scheme for 
computing steady supercritical flow about symmetrical air- 
foils and applied it to an ellipse for zero angle of attack. 
In this study, an algorithmic description of this new scheme 
is presented. Application to a symmetrical bicircular arc 
airfoil is also proposed. 

In Li and Holt's scheme, both Telenin's Method and the 
Method of Lines are used. The flow field before the shock 
is region 1. For transonic flow, singularity can be avoided 
by integrating the resulting ordinary differential equations 
away from the body. Region 2 contains the shock which will 
be located by shock fitting techniques. The shock divides 
region 2 into supersonic and subsonic regions and there is 
no singularity problem in this case. The Method of Lines 
is used in this region and it is advantageous to integrate 
the resulting ordinary differential equation along the body 
for shock fitting. 

Coaxial coordinates have to be used for the bicircular 
arc airfoil so that boundary values on the airfoil body can 
be taken with one direction of the coaxial coordinates fixed. 
To avoid taking boundary values at ±°° in the coaxial co- 
ordinary system, approximate analytical representation of 
the flow field near the tips of the airfoil is proposed. 
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1 . Introduction 

The major problem in transonic flow calculations is the 
nonlinearity in the governing equations and it causes the 
partial differential equations to change type within the 
solution domain, from elliptic in the subsonic region to 
hyperbolic in the supersonic region. The embedded super- 
sonic region ends with a shock, which causes an additional 
problem of handling the shock wave. 

There are two major categories of finite difference 
techniques developed to solve the transonic flow problem. 

One solves the transonic small disturbance equations and the 
other solves the full potential equations. 

In the first category, Murman and Cole (1971) developed 
the first efficient successive line overrelaxation (SLOR) 
method for the transonic small disturbance equations. They 
used a mixed finite difference system. The use of elliptic 
or hyperbolic difference formulas depends on whether the 
flow field is subsonic or supersonic. Shock capturing was 
used to locate the shock. Krupp and Murman (1972) extended 
the method to lifting airfoils and slender bodies. Later 
Ballhaus, Jameson and Albert (1978) developed an implicit 
approximate factorization algorithm for the solution of 
steady transonic small disturbance equations and found it 
to have a better rate of convergence than the SLOR algorithm. 

In the second category Steger and Lomax (1972) treated 
the transonic lifting airfoils by solving the full potential 
equations using SLOR. Jameson (1974) improved it by intro- 
ducing a rotated differencing scheme in which the direction 
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of upwind differencing is rotated to conform with the local 
flow direction. The most recent work in this category 
known to this author was done by Holst and Ballhaus, who 
solved the full potential equations in conservation form 
to ensure conservative shock capturing. They later used 
an arbitrary mesh for the method and obtained good results. 

Interpolation (semi-analytical) techniques that have 
been used in this problem include the Method of Integral 
Relations, Telenin's Method and the Method of Lines. 

However, the scheme developed by Li and Holt is the first 
one that can be generally applied to symmetrical airfoils 
and yet is able to locate the shock completely. 

In the Li and Holt scheme, the steady two-dimensional 
full potential equations are solved by both Telenin's 
Method and the Method of Lines. A doublet solution for 
flow past a closed body is used as the far field boundary 
condition. The strength of the doublet is a function of 
both the profile of the airfoil and the flow field. There- 
fore, it is an unknown and iterations on the strength of 
the doublet have to be done to obtain the correct boundary 
values on the airfoil. Jump conditions of the governing 
equations are applied across the shock wave so that it is 
perfectly sharp. Correct shock location is obtained by 
iterations and checks with the boundary conditions down- 
stream of the shock. The iterations mentioned above are 
two point boundary value problems and can be done efficiently 
by Powell's method. 
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Analytical representation of the flow field near the 
tips of the airfoil can be approximated by the subsonic 
small disturbance flow. Analytical representation, as 
described above, is required for the bicircular arc airfoil 
because one of the coordinates in the coaxial system approaches 
infinity at the two ends of the airfoil. These will be 
further elaborated upon in section 3.1. 
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2 . Formulation of the Problem 

In this study, a two-dimensional uniform flow past a 
symmetrical bicircular arc airfoil at zero angle of attack 
is considered. The embedded supersonic region near the 
maximum thickness section is as shown in Fig. 1. It is 
assumed that the flow is both irrotational and isentropic. 

In fact, it was shown by Li and Holt (1981) that for thin 
airfoils with subsonic free streams, the shock wave strength 
is sufficiently small and entropy changes can be neglected. 

2 . 1 Coaxial coordinates 

Before proceeding further into formulation, the coaxial 
coordinates system used in this study should be introduced. 
It will be easier to understand the coaxial coordinates 
system by referring to Fig. 2. In Milne -Thomson (1968), 
it is defined as 

z = ic cot h 5 , £ = 5 + in (2.1) 

£ and n, the coaxial coordinates are defined as 

c sinh n _ c sin £ ^ 2) 

x “ cosh n - cos £ ' * cosh n - cos £ 

2c is the chord length of the airfoil. When £ = constant, 

Eq. (2.2) is a circle whose center is the point (0, c cot £) 
with radius c cosec £. When n = constant, Eq . (2.2) is a 
circle whose center is the point (c coth n, 0) with radius 
c cosech n- The metric coefficient h is calculated to be 

-1 


h = c(cosh n - cos £) 


(2.3) 
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u and v are velocities in 5 and n directions, respectively. 


2 . 2 Equations of motion 


Continuity 


g|- (hpu) + (hpv) = 0 , 


(2.4) 


Irrotationality 


95 


(hv) - 


9n 


(hu) = 0 , 


(2.5) 


Bernoulli's equation 


H+i = g 2 - H o - *4x - 


( 2 . 6 ) 


H is the enthalpy per unit mass, q is the flow speed, and 
q max the maximum steady expansion speed. 

By assuming perfect gas with constant specific heat 

gives 


H = 


J CR. 


( Y— 1 ) P 

and isentropic flow which further gives 


(2.7) 


It follows that 


( 2 . 8 ) 


H _ p P o ^ ( P ] Y~1 

TT « ' ' 


H P P 
o ^o 


Therefore, Eq. (2.6) can be written as 

1 

V 1 


jl , 

w o H max 


(2.9) 


( 2 . 10 ) 
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If we express all variables in dimensionless form by dividing 
distances by c, velocities by *3 max » and the density by the 
stagnant density p q , retaining the same symbols for the non- 
dimensional variations, Eqs. 


(Puh) + (pvh) = 0 , 


9 (vh) J_ 

ac " an 


(uh) = 0 , 


,, 2 2 . y-1 

p = ( 1 - u - v ) 


( 2 . 11 ) 


( 2 . 12 ) 


(2.13) 


Equations (2.11), (2.12) and (2.13) are the three equations 

of motion with the three variables u,v and p. 
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3 . Numerical Methods 

Telenin's Method as applied to two-dimensional problems 
uses smooth interpolating polynomials to represent the 
unknowns in one of the independent variables. The governing 
partial differential equations are solved in their original 
form, and thus we avoid the algebra required in the Method 
of Integral Relations. The symmetry conditions of the present 
problem suggest the use of Fourier series of the form 

N 

u(£,n) = 2 a. (?) cos (i-1) n , (3.1) 

i=l 1 

N 

v(£,n) = £ a. U)sin(i-l)n / (3.2) 

i=l 

where N is the number of rays. In the present problem, the 
left most and right most rays do not fall on n = -°° and 
n =+oo d ue to the special handling of boundary conditions at 
these two regions. On jth ray 

N 

u. = u(£,n.) = £ a. ( 5) cos (i-l)n . . (3.3) 

3 ' J i -i 1 ' J 


The matrix {cos (i-1) n j } can be inverted to obtain a^, 

N 

a. = Z A . . u . , i=l,...,N , (3.4) 

i j =1 i] ] 

where {A i ^ } = {cos ( j-1) n i > 1 . 

From Eq. (3.1), 


(— ) 

v 3n ' i 


Substituting Eq. 


(— ) 


3u(£,n« ) n 

2 — = E a. (?) (1-i) sin (i-1) n 0 

3n ■ *. i * 


(3.4) into (3.5) yields 
N N 

£ ( £ A. .u.) (1-i) sin (i-1) r\ 
i=l j=i 1 3 3 *- 


(3.5) 


(3.6) 
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Equation (3.6) can also be written as 


3u. 


N N 


N 


<$?>* '•! A i j(l-i)si n (i-l)n l )u j -Z P tjUj (3.7) 

] J. i x 3 


Similarly, 


3v x 


N 


whe re 


and 


( a nh “ j f 1 G aj v j ' 1 = 1 ' 2/ ***' N ' 


N 

G 0 . = £ B. . (i-l)cos (i-l)n„ 

X'J x=r]_ ^ 


(3.8) 


(3.9) 


{B^} = {sin( j-l)r)j} _1 . 

Equations (3.7) and (3.8) are used as interpolations of u and 
v in the n direction; therefore, the resulting ordinary dif- 
ferential equations can be integrated in the £ direction. 

The Method of Lines is very similar. Instead of interpolating 
polynomials, the n derivatives are approximated by three- 
point or five-point finite difference schemes. 

We can obtain the expressions for — and rr from Eqs. 

dt, a c, 

(2.11), (2.12) and (2.13). From Eq. (2.12) 


v, 3v j. ,, - v, 9u j. „ 3h 

h 3? +V H~ h H +U 3W • 


(3.10) 


Therefore, 


5v _ 1 r. 3u 3h 3h 1 

3f = ht h ^ + U 37f- V 3S ] 


3n 


But 


3h 

3n 


-r— = — n sinh n 


3h _ ,2 . 

— - -h sin ? 


(3.11) 


(3.12) 


and gives 


^ - uh sinh n + vh sin £ 

3£ 3n 


(3.13) 
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Similarly, from Eqs. (2.11) and (2.13), 

3u _ P 
3C " Q 

where 

P = 2uv 4~: - uh sinh n +vh sin £ 


2 2 

+ (y-l) (1-u -v ) [uh sin£+vh sinh n 


, 0 3u . 3v 

+ 2v u aTT + v 377 ' 


Q = (Y-l) (l-u 2 -v 2 ) - 2u 2 


(3.15) 



(3.15a) 

(3.15b) 


Therefore, the singular ellipse is obtained from Q=0, or 


where 



2 = Xl_ 

Y+l * 


(3.16) 

(3.17) 


All points on the ellipse lie outside the sonic circle, 

* 

except for v=0, u= q . Therefore, it is apparent that for 
the supercritical flow as in the present problem, no 
singularity will be encountered when integrating in the 
£ direction, or away from the body. Therefore, for region 1 
(see Fig. 4), we must integrate the resulting ODE's from 
the body. 

The expressions for and can be similarly obtained 
when we interpolate in the £ direction and integrate in the 
H direction. This is applied to region 2 where supersonic 
and subsonic regions are separated by the shock. In this 
case, the singular ellipse is 


u 2 + 


2 

V 

*2 


1 


(3.18) 
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Any critical or supercritical flow has a point on the body 

* 

with u=0 and v =q , which is a point on the ellipse given 
by (3.18). That shows that integration in the n direction 
in region 1 always leads to at least one singularity at 
sonic points on or near the body. 

3 . 1 Boundary conditions 

For the solution represented by Eqs. (2.11), (2.12) 
and (2.13), we need the boundary values for u and v on the 
body and at n = - 00 . But to fit the shock, we also need the 
boundary conditions at n = + 00 so that iterations can correct 
the position of the shock by checking with values of u and 
v at n = +«> . 

On the body, the normal velocity is zero, or 

u = 0 for £ = C Q . (3.19) 

Since the tangential component v is not known, the far field 
boundary conditions are needed. Murman and Cole (1971) 
derived an analytical solution for the far field by using 
the transonic small-distixrbance equation 

[K<f> x -Js(Y+l)^] x + <f> yy = 0 (3.20) 

with the variables and parameters defined by 

y = <S 1/3 y , K = ( 1 - M^) / 6 2/3 . (3.21) 

The far field they obtained is that of the usual doublet 
for a closed body 
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<J> (x,y) 


D x 

27TK* 5 (x 2 + Ky 2 ) 


(3 


where the doublet strength is 


D = 2 


-1 


F (x) dx +H (Y+l ) 


{<3^(x,y)} dxdy . (3 


The perturbation velocities 



(3 


There fore 


with 


q ' 

x 


2 ~ 2 
D (-x +Ky ) 

V — 5 Zo 3 

2ttK 2 (x +Ky^) 




_ x y 

(x 2 + Ky 2 ) 2 


q x 

— = 1 
U 

oo 


6 2 / 3 q' 

n x 


i. 

u 


6q ' 


(3 


The flow velocities expressed in the £ and n directions 
are given by 


u 



_ _1 r ~ 1 

h L sinh n sin C q x 


cos E, cosh n 


osh n - 1 ^ ^ 


v = ^ (hn) = £ [ 


h L 1 - cosh n cos E, 


sin £ sinh n H y 


q,J ' 


(3. 


The coordinate n goes to ±°° at the two ends of the 


. 22 ) 


.23) 


.24) 


(3.25) 


.26) 


27) 


28) 


airfoil, but remains finite except at a small distance 
very close to the tips. Analytical representation of the 
flow field in these two regions has to be found so that 
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boundary conditions can be taken at finite values of n* 
Figure 3 shows the velocity distribution on the ellipse 

with thickness ratio 6 =0.4 for M =0.65 (Li and Holt 1981). 

00 

The velocity distribution near the two ends appears to be 
symmetrical. Moreover, the velocity in this region is of 
Mach number less than 0.8. Therefore, subsonic linearized 
potential equation (Shapiro 1953) can be applied and used 
as boundary conditions. The basic equation is 

(1 -M 2 ) + li| = o . (3.29) 

Ox Oy 

It can also be written as 

= 0 . (3.30) 

Ox' Oy 

where 


x' = - . (3.31) 

a -«» 

If we represent the airfoil profile as 

y = fif(x') , ~ C - ' < x 1 < c — (3.32) 

VI -M 2 /I -M 2 

CO oo 

and the potential as 


<p = Ux' + 6 4> 1 , (3.33) 

where 6 is the thickness ratio of the airfoil, then the 
subsonic velocity distribution is known when <j>^ is 
determined. 

We can use plane source with strength r ( x ) so as to 
satisfy the boundary conditions of the subsonic small- 
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disturbance flow. Hence 


♦l = 


£n/(x' - X ) 2 + y 2 d X 


where 


c' = 


/I -M 


On the surface of the thin airfoil 


!£= «f(x') 


If we expand ||- and about y =0, substituting them into 

(3.36) and compare the first order (in 6 ) terms, we get 


(x- ,0) = u f ' (x') . 


From Eqs. (3.34) and (3.37), and letting t = x — 

Xj ™ x ^ 

= i r U'+yt) „„ 


)y 2 tt 


1+t ' 

y 


Letting y 0, 


5 ^r r (x'+yt) 

-CO 217 l + t‘ 


- H r(x') . 


Therefore, we finally obtain 


T(x') = 2u f* ( x 1 ) 


and <f>^ can thus be determined, 


_ 8 <p 


_ 9 ({> 


‘x 3x ' H y 9y 


(3.42) 
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are then used as boundary conditions for the region where 
n goes to ±°°. 

3 . 2 Jump conditions 

The shock wave on the downstream side of the supersonic 
flow over the maximum thickness region of the body is handled 
by shock- fitting in the present problem. The shock wave is 
modelled by a jump discontinuity in the solution and the 
jump conditions are satisfied exactly so that the shock 
wave is perfectly sharp. Instead of the usual Rankine- 
Hugoniot relations, the jump conditions can be derived from 
the equations of motion in the present problem. Applying 
the two-dimensional form of the divergence theorem to 
Eqs. (2.4) and (-2.5), we get 

<puh> (dn) - <pvh> (d£) = 0 (3.43) 

and 

<vh> (dn) + <uh> (d£) = 0 (3.44) 

s s 

where < > denote a jump in the quantity across the shock and 
subscript s denotes an element in the shock surface. Equa- 
tions (3.43) and (3.44) can also be written as 


and 


where 


<puh>n ' - <pvh> = 0 
s 

<vh>rig + <uh> = 0 , 

n; = (^j g , the shock wave angle 


(3.45) 

(3.46) 


(3.47) 


The metric coefficient h is continuous and has the same 
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value on both sides of the shock surface and can be dropped 
from Eqs. (3.45) and (3.46). Hence the final form of the 
jump conditions is 

<P u>n g - <p v> = 0 , (3.47) 

<v>Hg + <u> = 0 . (3.48) 

3. 3 Implementation of the numerical scheme 

Telenin's Method and the Method of Lines as applied to 
elliptic partial differential equations solve a Dirichlet 
problem as a Cauchy problem. It is inherently unstable with 
respect to the prescribed data. This phenomenon is known 
as Hadamard instability. Jones, South and Klunker (1972) 
encountered Hadamard instability in applying the Method of 
Lines and found growth in error proportional to exp(N£), 
where N is the number of rays and £ the direction of 
integration. To have sufficient number of rays to repre- 
sent the variables near the body, we are thus restricted 
in the direction of integration, £ . Following Li and 
Holt (1981) , the present problem can be solved in two stages 
to overcome the difficulty. In the first stage, a very 
coarse representation of the variables is used, which 
enables us to integrate the equations away from the body 
to the far field without instability problems. We will 
obtain a supercritical shock free flow, which is unstable 
and not likely to occur in practical situations. The coarse 
solution provides a fairly good representation of the flow 
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field away from the body where the flow is smooth. The 
coarse solution at an intermediate value of £, say £^, 
is used as the outer boundary condition. Therefore, a 
larger number of rays can be used for the refined solution 
near the body. We have to choose £^ so as to avoid having 
rays of constant £ with values close to £^ to pass through 
the sonic line that will cause difficulty when integrate 
through this line. The two regions of integration are 
shown in Fig. 4. 

An algorithmic description of the numerical scheme 
proposed for the present problem is presented as follows. 
First, we need to 

(a) obtain the simultaneous ordinary differential equa- 
tions (Eqs. (3.13), (3.15)) by interpolating the 

variables in one direction. 

(b) estimate the tangential velocity on the surface of 
the airfoil. 

(c) obtain an analytical representation of the flow field 
near n =±°° as boundary conditions. 

(d) estimate the doublet strength D for the far field 
doublet solution (Eq. (3.22)). 

The computational procedures are then 

1. integrate the ode's with large steps in £ and n . 

Use Telenin's Method for the coarse solution to 
obtain a coarse solution. Integration in the second 
quadrant only is required as the shock free flow is 
symmetrical. 
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2 . integrate in £ direction and compare with the far 
field doublet solution. 

(a) correct the estimates for tangential velocities and 
repeat until the difference between the far field 
solutions falls within tolerance. It can be done 
effectively with Powell's Method which will be 
briefly explained in section 3.3a. 

(b) recalculate D with the results obtained in 2(a) 
and repeat 2 and 2(a) until values of D converge. 

3. the coarse solution hence obtained at is used as 
boundary conditions for the refined solution. 

4. divide the flow field into two regions (Fig. 4) and 
refine the mesh. 

5. repeat the procedure 2 for the refined solution in 
region 1. 

6. For region 2, the ode's are integrated in the n 
direction. 

7. estimate the location of the shock. The shock wave 

is normal on the surface. The jump conditions ((3.47) , 
(3.48)) are fully specified when the shock location is 
known . 

8. integrate in the n direction, apply jump conditions 
across the shock. Check with the downstream boundary 
conditions . 

9. correct the shock location with Powell’s Method. 

10. Lastly, recalculate doublet strength D with the solu- 
tion obtained and repeat the whole procedure . The 
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whole computation is complete when values of D at 
successive iterations agree to within the prescribed 
tolerance . 

The numerical scheme is briefly summarized by the flow 
chart in Fig. 5. 


3 . 3a Powell's Method 

Powell's Method is essentially that of least squares 
minimization. A concise description of Powell's method can 
be found in section 6.9 of Numerical Methods in Fluid 
Dynamics by Holt [1984]." This method is best explained 
by example. 

In the procedure 2 of the numerical procedure, the 
difference between the two far field solutions depends 
on the initial guess of the tangential velocity, Vj . 

Then the method minimizes 


N 2 


(3.49) 


with respect to V^, . . . , V N - is minimized by making 

changes to Vj according to the direction 6V given by 


N N 
Z { Z 


w 1 = - z e k ~ 

j=l k=l 3V i 8V j 2 3 k=l k 3V i 


3e, 


i = 1 , 2, . . . ,N 

(3.50) 


New values of V are given by 


V = V , , + X6V , (3.51) 

old 

2 

in which X is chosen such that is minimized along the 
direction <5V. The required X can be chosen by evaluating 
at different values of X. 
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4 . Conclusion 

The application of the composite numerical scheme 
developed by Li and Holt (1981) to a bicircular arc air- 
foil is proposed. As coaxial coordinates have to be used, 
analytical representation of the flow field near the two 
tips of the airfoil is required as boundary conditions. 

The analytical representation can be easily constructed 
by assuming a subsonic small-disturbance flow at the two 
ends of the airfoil. An algorithmic description of the 
numerical scheme is also presented. 


February 1989 
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Fig. 1(a) Unstable Supercritical Shock Free Flow Field. 
Upper Half Plane. 



Fig. 1(b) Typical Supercritical Flow Field. Upper 
Half Plane. 



Fig. 2 The construction of coaxial coordinates, n =^n(r 2 /r 1 ), 
£ 0 =0^ “02* The ha l f “Pl ane with negative y is used. 
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to initiate with 
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to compute coarse solution 


integrate ODE's in E, direction 

iterate for correct V^'s on surface with Powell's method 
recalculate D and repeat until it converges to D 


coarse solution at used 
as boundary conditions 


to compute refined solution 

integrate ODE's in £ direction for region 1 
iterate for correct Vj_'s on surface with Powell's method 
integrate ODE ' s in direction for region 2 
apply jump conditions across the shock 
shock location is adjusted with Powell's method and 
checked with the downstream BC's 


recalculate D 
obtain D 


Is (D^-Djl < tolerance? 


D i = D f 



Fig. 5 Flow chart of the numerical scheme 






